<!DOCTYPE html>
<html>

<head>
  <title>使用leaflet画椭圆，弧线，根据起点经纬度，速度，时间，方向，计算终点经纬度-方法一</title>
  <meta charset="utf-8" />
  <meta name="viewport" content="width=device-width, initial-scale=1.0">
  <link rel="stylesheet" href="https://unpkg.com/leaflet@1.6.0/dist/leaflet.css"
    integrity="sha512-xwE/Az9zrjBIphAcBb3F6JVqxf46+CDLwfLMHloNu6KEQCAWi6HcDUbeOfBIptF7tcCzusKFjFw2yuvEpDL9wQ=="
    crossorigin="" />
  <script src="https://unpkg.com/leaflet@1.6.0/dist/leaflet.js"
    integrity="sha512-gZwIG9x3wUXg2hdXF6+rVkLF/0Vi9U8D2Ntg4Ga5I5BZpVkVxlJWbSQtXPSiUTtC0TjtGOmxa1AJPuV0CPthew=="
    crossorigin=""></script>
</head>
<style>
  html,
  body {
    height: 100%;
  }

  #mapid {
    height: 100%;
    overflow: hidden;
  }
</style>

<body>
  <div id="mapid"></div>
  <script>
    /**
     * 算法背景：公司需要一个模型，求污染扩散的面积，给出起始的经纬度坐标，风速，风向，时间，根据这四个参数画出一个椭圆扩散的图形，起点到中间位置为重污染红色，中间位置到3/4位置为中污染黄色，3/4位置到终点为轻污染绿色。
     * destinationVincenty()方法计算终点位置
     * 依次画出起点到终点，起点到3/4点，起点到中点的线段
     * calcCoorArr()方法画出起点到各点的弧线，一共六条
     */
    var map = L.map('mapid').setView([37.043797, 114.44556], 12);
    L.tileLayer('http://wprd0{s}.is.autonavi.com/appmaptile?x={x}&y={y}&z={z}&lang=zh_cn&size=1&scl=1&style=7', {
      attribution: '&copy; 高德地图',
      maxZoom: 24,
      minZoom: 4,
      subdomains: "1234"
    }).addTo(map);

    drawEllipse({ lon: 114.44556, lat: 37.043797 }, 4, 1, 270)

    /**
    * 给定位置经纬度方向时间速度画扩散范围
    * @lonlat  {Object}  经纬度  {lon:0,lat:0}
    * @speed  {number} 速度 m/s
    * @time  {number} 时间 h
    * @direction {number}   方向
    */
    function drawEllipse(lonlat, speed, time, direction) {
      //定义起点
      var lonlat = lonlat
      //定义划线起点red,yellow,green
      var latlngsR = [[lonlat.lat, lonlat.lon]], latlngsY = [[lonlat.lat, lonlat.lon]], latlngsG = [[lonlat.lat, lonlat.lon]]
      var speed = speed;  // 速度 m/s
      var time = time;   // 时间 hour
      var distence = speed * time * 3600;
      var direction = direction - 180 //角度
      var VincentyConstants = {
        a: 6378137,
        b: 6356752.3142,
        f: 1 / 298.257223563
      }
      latlngsR.push(destinationVincenty(lonlat, direction, distence / 2))
      latlngsY.push(destinationVincenty(lonlat, direction, distence * 3 / 4))
      latlngsG.push(destinationVincenty(lonlat, direction, distence))
      //使用calcCoorArr方法计算出路径点
      var latlngsRB = calcCoorArr(latlngsR[0][1].toString() + ',' + latlngsR[0][0].toString(), latlngsR[1][1].toString() + ',' + latlngsR[1][0].toString(), 100, 3);
      var latlngsRE = calcCoorArr(latlngsR[1][1].toString() + ',' + latlngsR[1][0].toString(), latlngsR[0][1].toString() + ',' + latlngsR[0][0].toString(), 100, 3);
      var latlngsYB = calcCoorArr(latlngsY[0][1].toString() + ',' + latlngsY[0][0].toString(), latlngsY[1][1].toString() + ',' + latlngsY[1][0].toString(), 100, 3);
      var latlngsYE = calcCoorArr(latlngsY[1][1].toString() + ',' + latlngsY[1][0].toString(), latlngsY[0][1].toString() + ',' + latlngsY[0][0].toString(), 100, 3);
      var latlngsGB = calcCoorArr(latlngsG[0][1].toString() + ',' + latlngsG[0][0].toString(), latlngsG[1][1].toString() + ',' + latlngsG[1][0].toString(), 100, 3);
      var latlngsGE = calcCoorArr(latlngsG[1][1].toString() + ',' + latlngsG[1][0].toString(), latlngsG[0][1].toString() + ',' + latlngsG[0][0].toString(), 100, 3);
      //绿色面
      L.polygon([...latlngsGB, ...latlngsGE], { color: 'green', fillOpacity: 0.8, fill: true, fillColor: 'green' }).addTo(map);
      //黄色面
      L.polygon([...latlngsYB, ...latlngsYE], { color: 'yellow', fillOpacity: 0.6, fill: true, fillColor: 'yellow' }).addTo(map);
      //红色面
      L.polygon([...latlngsRB, ...latlngsRE], { color: 'red', fillOpacity: 0.4, fill: true, fillColor: 'red' }).addTo(map);

      /**
      *Calculate destination point given start point lat/long (numeric degrees),
      * bearing (numeric degrees) & distance (in m).
      * 给定位置经纬度方向距离计算目的地点经纬度
      * @lonlat  {Object}  经纬度  {lon:0,lat:0}
      * @brng  {number} 度数 0—360 0：北方  180：南方
      * @dist  {number} 距离(m)
      * @return {Array}   经纬度
      */
      function destinationVincenty(lonlat, brng, dist) {
        var ct = VincentyConstants;
        var a = ct.a, b = ct.b, f = ct.f;
        var lon1 = lonlat.lon * 1;  //乘一（*1）是为了确保经纬度的数据类型为number
        var lat1 = lonlat.lat * 1;
        var s = dist;
        var alpha1 = rad(brng);
        var sinAlpha1 = Math.sin(alpha1);
        var cosAlpha1 = Math.cos(alpha1);
        var tanU1 = (1 - f) * Math.tan(rad(lat1));
        var cosU1 = 1 / Math.sqrt((1 + tanU1 * tanU1)), sinU1 = tanU1 * cosU1;
        var sigma1 = Math.atan2(tanU1, cosAlpha1);
        var sinAlpha = cosU1 * sinAlpha1;
        var cosSqAlpha = 1 - sinAlpha * sinAlpha;
        var uSq = cosSqAlpha * (a * a - b * b) / (b * b);
        var A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
        var B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));

        var sigma = s / (b * A), sigmaP = 2 * Math.PI;
        while (Math.abs(sigma - sigmaP) > 1e-12) {
          var cos2SigmaM = Math.cos(2 * sigma1 + sigma);
          var sinSigma = Math.sin(sigma);
          var cosSigma = Math.cos(sigma);
          var deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
            B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
          sigmaP = sigma;
          sigma = s / (b * A) + deltaSigma;
        }

        var tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1;
        var lat2 = Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,
          (1 - f) * Math.sqrt(sinAlpha * sinAlpha + tmp * tmp));
        var lambda = Math.atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1);
        var C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
        var L = lambda - (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));

        var revAz = Math.atan2(sinAlpha, -tmp);  // final bearing
        var lon_destina = lon1 * 1 + deg(L);
        var lonlat_destination = [deg(lat2), lon_destina];

        return lonlat_destination;
      }

      /**
      * 度换成弧度
      * @param  {Float} d  度
      * @return {[Float}   弧度
      */
      function rad(d) {
        return d * Math.PI / 180.0;
      }

      /**
       * 弧度换成度
       * @param  {Float} x 弧度
       * @return {Float}   度
       */
      function deg(x) {
        return x * 180 / Math.PI;
      }
      /**
       * 画曲线
       * @point_start  {String} 起点坐标 "114.30911,30.600052"
       * @point_end {String} 起点坐标  "114.30911,30.600052"
       * @num {Number} 坐标点个数
       * @w {Number} 弯曲度
       * @return {Array} 曲线坐标点
       */
      function calcCoorArr(point_start, point_end, num, w) {
        //第一步把坐标字符串转为对象,为了方便计算转为了数字格式
        var p_start = { x: parseFloat(point_start.split(",")[0]), y: parseFloat(point_start.split(",")[1]) };
        var p_end = { x: parseFloat(point_end.split(",")[0]), y: parseFloat(point_end.split(",")[1]) };
        //此处敲黑板，是任务的第二大难点，求出相对的第三个点，用于固定曲线的弯曲度，下面公式是已知三角形两点坐标和两个坐标点的夹角求第三点坐标，两个夹角我们是自定义任意值，不过不要加起来超过180度
        let x3 = (p_start.x * w + p_end.x * w - p_start.y + p_end.y) / (2 * w)
        let y3 = (p_start.y * w + p_end.y * w - p_end.x + p_start.x) / (2 * w)
        var p_crt1 = { x: x3, y: y3 };
        var p_crt2 = { x: x3, y: y3 };
        let paths = [];
        for (let i = 0; i < num + 1; i++) {
          let t = i / num;
          var _matrix1 = [1, t, t * t, t * t * t];
          var _matrix2 = [
            [1, 0, 0, 0]
            , [-3, 3, 0, 0]
            , [3, -6, 3, 0]
            , [-1, 3, -3, 1]
          ];
          var _matrix3 = [
            [p_start.x, p_start.y]
            , [p_crt1.x, p_crt1.y]
            , [p_crt2.x, p_crt2.y]
            , [p_end.x, p_end.y]
          ];
          var _matrix_tmp = [
            _matrix1[0] * _matrix2[0][0] + _matrix1[1] * _matrix2[1][0] + _matrix1[2] * _matrix2[2][0] + _matrix1[3] * _matrix2[3][0]
            , _matrix1[0] * _matrix2[0][1] + _matrix1[1] * _matrix2[1][1] + _matrix1[2] * _matrix2[2][1] + _matrix1[3] * _matrix2[3][1]
            , _matrix1[0] * _matrix2[0][2] + _matrix1[1] * _matrix2[1][2] + _matrix1[2] * _matrix2[2][2] + _matrix1[3] * _matrix2[3][2]
            , _matrix1[0] * _matrix2[0][3] + _matrix1[1] * _matrix2[1][3] + _matrix1[2] * _matrix2[2][3] + _matrix1[3] * _matrix2[3][3]
          ];
          var _matrix_final = [
            _matrix_tmp[0] * _matrix3[0][0] + _matrix_tmp[1] * _matrix3[1][0] + _matrix_tmp[2] * _matrix3[2][0] + _matrix_tmp[3] * _matrix3[3][0]
            , _matrix_tmp[0] * _matrix3[0][1] + _matrix_tmp[1] * _matrix3[1][1] + _matrix_tmp[2] * _matrix3[2][1] + _matrix_tmp[3] * _matrix3[3][1]
          ];

          var _res_point = [_matrix_final[1], _matrix_final[0]];
          paths.push(_res_point);
        }
        return paths;
      }

    }
  </script>
</body>

</html>